// -----------------------------------------------------------------------------
// M A R I T I M E   R E S E A R C H   I N S T I T U T E   N E T H E R L A N D S
// -----------------------------------------------------------------------------
// Copyright (C) October 16, 2009 MARIN - All rights reserved - http://www.marin.nl
// -----------------------------------------------------------------------------
// Program  : various
// Module   : ???
// File     : Qshep2d.cpp
// -----------------------------------------------------------------------------

#include <new>
#include <math.h>
#include "Qshep2d.h"
#include "qshep2d_c.h"
#include "assertplus.h"

// -----------------------------------------------------------------------------
// function Qshep2d::Qshep2d
// -----------------------------------------------------------------------------
/// constructor
/** detailed description */
// -----------------------------------------------------------------------------

Qshep2d::Qshep2d()
: m_n(0)
, m_x(0)
, m_y(0)
, m_z(0)
, m_lcell(0)
, m_lnext(0)
, m_xmin(0)
, m_ymin(0)
, m_dx(0)
, m_dy(0)
, m_rmax(0)
, m_rsq(0)
, m_a(0)
, m_nr(0)
{
}


// -----------------------------------------------------------------------------
// function Qshep2d::~Qshep2d()
// -----------------------------------------------------------------------------
/// destructor
/** detailed description */
// -----------------------------------------------------------------------------

Qshep2d::~Qshep2d()
{
    delete [] m_x;
    delete [] m_y;
    delete [] m_z;
    delete [] m_lcell;
    delete [] m_lnext;
    delete [] m_rsq;
    delete [] m_a;
}


// -----------------------------------------------------------------------------
// function Qshep2d::setKnownValues
// -----------------------------------------------------------------------------
/// enters known xyz coordinates into interpolation object
/// @return success
// -----------------------------------------------------------------------------

bool Qshep2d::setKnownValues(const Vertices& xy, const Values& z)
{
    // fill data structures for qshep2
    m_n = (int)(xy.size());


    int nq = 13; // the number of data points to be used in the least
                 // squares fit for coefficients defining the nodal functions Q(K).
                 // A highly recommended value is NQ = 13.
                 // NQ must be at least 5, and no greater than the minimum of 40 and N-1.
    int nw = 19; // the number of nodes within (and defining) the radii
                 // of influence R(K) which enter into the weights W(K).  For N
                 // sufficiently large, a recommended value is NW = 19.
                 // NW must be at least 1, and no greater than the minimum of 40 and N-1.
    m_nr = 3;    // the number of rows and columns in the cell grid
                 // defined in subroutine STORE2.  A rectangle containing the nodes
                 // is partitioned into cells in order to increase search efficiency.
                 // NR = SQRT(N/3) is recommended.  NR must be at least 1.

    ASSERT(m_n >= 6, "use at least 6 data points");
    int maxnqw = std::min(40, m_n - 1);
    if (nq > maxnqw)
    {
        nq = maxnqw;
    }
    if (nw > maxnqw)
    {
        nw = maxnqw;
    }
    m_nr = int(sqrt(m_n / 3.0));

    delete m_x;
    m_x = new double[m_n];
    delete m_y;
    m_y = new double[m_n];
    delete m_z;
    m_z = new double[m_n];
    delete m_lcell;
    m_lcell = new int[m_nr * m_nr];
    delete m_lnext;
    m_lnext = new int[m_n];
    delete m_rsq;
    m_rsq = new double[m_n];
    delete m_a;
    m_a = new double[5 * m_n];
    for (int i = 0; i < m_n; ++i)
    {
        m_x[i] = xy[i].x;
        m_y[i] = xy[i].y;
        m_z[i] = z[i];
    }

    //  Call QSHEP2 to define the interpolant Q to the data
    int ier = 0;
    qshep2(m_n, m_x, m_y, m_z, nq, nw, m_nr, m_lcell, m_lnext,
           &m_xmin, &m_ymin, &m_dx, &m_dy, &m_rmax, m_rsq, m_a, &ier);
    if (ier != 0)
    {
        printf("%s %8d\n", "  Error in QSHEP2D, IER = ", ier);
        return false;
    }
    return true;
}


// -----------------------------------------------------------------------------
// function Qshep2d::getUnknownValue
// -----------------------------------------------------------------------------
/// fills in z-value at given xy coordinate
// -----------------------------------------------------------------------------

double Qshep2d::getUnknownValue(const Vertex& xy)
{
    double px = xy.x;
    double py = xy.y;
    double pz = qs2val(px, py, m_n, m_x, m_y, m_z, m_nr, m_lcell, m_lnext,
                       m_xmin, m_ymin, m_dx, m_dy, m_rmax, m_rsq, m_a);
    return pz;
}


// -----------------------------------------------------------------------------
// function Qshep2d::getUnknownValues
// -----------------------------------------------------------------------------
/// fills in z-values at given xy coordinates
// -----------------------------------------------------------------------------

bool Qshep2d::getUnknownValue(const Vertex& xy, double *pZ, double *pDX, double *pDY)
{
    double px = xy.x;
    double py = xy.y;
    int ier;
    qs2grd(px, py, m_n, m_x, m_y, m_z, m_nr, m_lcell, m_lnext,
           m_xmin, m_ymin, m_dx, m_dy, m_rmax, m_rsq, m_a, pZ, pDX, pDY, &ier);
    if (ier != 0)
    {
        printf("%s %6d\n", "  Error in QS2GRD, IER = ", ier);
        return false;
    }
    return true;
}
